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Abstract 

We show that Kolmogorov multipliers in turbulence cannot be sta- 
tistically independent of others at adjacent scales (or even a finite 
range apart) by numerical simulation of a shell model and by theory. 
As the simplest generalization of independent distributions, we sup- 
pose that the steady-state statistics of multipliers in the shell model 
are given by a translation-invariant Gibbs measure with a short-range 
potential, when expressed in terms of suitable "spin" variables: real- 
valued spins that are logarithms of multipliers and XY-spins defined 
by local dynamical phases. Numerical evidence is presented in favor 
of the hypothesis for the shell model, in particular novel scaling laws 
and derivative relations predicted by the existence of a thermodynamic 
limit. The Gibbs measure appears to be in a high-temperature, unique- 
phase regime with "paramagnetic" spin order. 
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1 Introduction 



In a famous paper on the local structure of turbulence of incompressible 
fluids [1], A. N. Kolmogorov in 1962 considered inertial-range multipli- 
ers denned by ratios of velocity increments Wij (£,£') = Si^Vj/Si g/Vj. Here 
5i iVj(x) = Vj(x + fej) — Vj{pc) is the increment of the jth component of the 
velocity vector v along the unit vector e^. Kolmogorov hypothesized that, 
at very high Reynolds number, these multipliers should have distributions 
which are universal functions only of the scale-ratio £/£' and not of the ab- 
solute scale. He postulated further that multipliers corresponding to widely 
separated scales should be statistically independent. See also [2]. R. Benzi, 
L. Biferale, and G. Parisi [3] have made a more precise and quite remarkable 
hypothesis in the context of a "shell model" of turbulence. The latter are 
quadratically nonlinear dynamical systems for variables in a finite number 
N of shells with wavenumbers k n = X n ko, n = 1,...,N with A > 1. For 
example, the SABRA model [4] has one complex mode u n per shell which 
obeys a dynamical equation 

i(ak n+ iu n+2 u* n+l + bk n u n +\Un-\ 
-cfc n _iu n _iu n _ 2 ) - vk 2 n u n + /„. (1) 

([3] considered a slightly different model.) Here * stands for complex conju- 
gation, f n is a forcing term which is restricted to the first few shells and v 
is the "viscosity". When a + b + c = 0, equation (1) satisfies conservation of 
"energy" E = \ J2 n \u n \ 2 and "helicity" H = £ n (a/c) n |« n | 2 in the v -» 
limit, analogous to the quadratic invariants of the inviscid Euler equations. 
Then, for large iV and high Reynolds number Re = \J (\ui\ 2 ) / (vko), the 
authors of [3] hypothesized that the variables u n , n = 1, N should have 
a steady-state statistical distribution given by a Gibbs measure 

P(u!,...,u N ) oc exp[-^$ n (« n ,u n _i,u n+ i, ...)]. (2) 

n 

Because turbulence is a dissipative state, far from thermodynamic equilib- 
rium, the potentials <I> n have nothing to do with the inviscid invariants, even 
as v — > 0. For shellnumbers n in the long inertial range, 1 < n C JV, [3] 
supposed that the potential $ n becomes a universal function independent 
of n. They supposed further that the potential $ is a sufficiently short-range 
function of "the ratios between the n's and their angles" [3]. The authors 
of [3] used an "infinite-temperature" model with independent multipliers to 
predict the scaling exponents. It is usually assumed that such an approxima- 
tion is qualitatively correct. However, we show that the multipliers cannot 
be strictly independent and that the correlations are essential. We also give 
evidence for the Gibbs hypothesis, based on theoretical analysis and direct 
numerical simulation of the shell-model dynamics (1). 



du n 
dt 
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2 Theoretical Considerations 



First, let us give a more precise form to the hypothesis. In the shell model, 
we introduce an amplitude p n and a phase 9 n for each shellnumber n, via 
the polar decomposition u n = k n 1 ^ ■ p n e ten . Following [3], we separate out 

— 1/3 

the Kolmogorov 1941 scaling factor k n . If w n = p n / p n -i is the multiplier 
defined by [3] and A n = — 9 n + Q n -\ + &n-2 is the dynamical phase factor 
[3, 4], then we define 

a n = ln(w n ), U n = exp(zA n ) (3) 

where a n is a local "slope" and U n is an "XY-spin" or 2-dimensional rotator 
spin. The definition of a n is motivated by the observation that In p n = 
J2k=i a k is then a "total spin" . Alternatively, h n = In p n can be viewed as a 
"height function" , as in equilibrium models of surface roughness. Then, our 
hypothesis is that the distribution of these "spin variables" £ n = (a n , U n ) is a 
translation-invariant Gibbs measure with a short-range potential <!>, ignoring 
finite-size effects from the forcing and dissipation ranges of shellnumber n. 
The potential is expected to be, at least, absolutely summable ([5], section 
2.1), which guarantees its uniqueness up to physical equivalence. However, 
the numerical evidence presented below suggests that the interactions are 
not merely summable, but indeed quite rapidly decaying. 

On the other hand, the potentials cannot have a strictly finite range 
R, i.e. vanishing for any set of spin variables containing pairs £ n ,£ n ' with 
\n — n'\ > R. In particular, the assumption of zero-range interactions, 
R = 0, or independent spins which was made made by [3] is ruled out. Us- 
ing exact constraints from the dynamics, we show in Appendix 1 that the 
assumption of independent spins leads to the deterministic K41 fixed-point 
u n ~ —i(e/k n ) 1 ^ 3 as the only statistically stable solution. In reality this 
solution is dynamically unstable [7]. The false assumption of independence 
stabilizes this solution and prevents intermittency corrections from devel- 
oping. This is plausible, since intermittency is known to arise in the shell 
models from "burst" solutions which exhibit long-range coherence in both 
a n and U n over many shells (e.g. see Fig. 4 in [8]). Furthermore, the K41 
fixed-point can be shown to be stabilized by assuming any finite-range po- 
tential between spins [6], so that the stationary measure cannot be Gibbs 
with any potential of strictly finite-range. 

Another argument against independence of spins is that this would imply 
the quadratic equation 

a/4 + bp 3 + c = (4) 

for the multiplier moment p p = (w p ), p = 3. Cf. [3], Eq.(18). This is anal- 
ogous to the "4/5-law" of fluid turbulence . Note that, in the independent 
spin approximation, the structure-function scaling exponents are given by 
Cp = (p/3) — logA Pp- There are two roots of the quadratic equation, ^3 = 1 
and P3 = c/a. The first solution gives constant mean energy flux but zero 
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helicity flux, while the second gives the opposite. Thus, no joint cascade 
of energy and helicity is possible in the independent spin approximation, 
contrary to observations [9]. Furthermore, the second solution violates the 
realizability inequality [13 > 0, when c/a < and the second invariant is 
truly "helicity-like" . These conclusions do not remain true for a potential of 
finite but non-zero range. The quadratic equation still holds for a nearest- 
neighbor potential or Markov chain approximation, where now ^3 is the 
principal eigenvalue of a "transfer matrix" [6]. As shown in Appendix 2 of 
the present paper, joint cascades of energy and helicity are permitted in a 
Markov chain model, if /i.3 = 1 and the subleading eigenvalue /U 3 = c/a. Fur- 
thermore, a concrete Markov chain model is constructed in the appendix to 
show that this situation may be realized. While the Gibbs measure cannot 
be exactly nearest-neighbor, we believe that this may be a good working 
approximation. 

3 Numerical Results 

We now present our simulation results for the SABRA model, with standard 
choice of parameters A = 2, ko = 2~ 4 , a = 1, b = c = —1/2. We performed 
two sets of simulations, one with N = 22, u = 1CT 7 and the second with 
N = 26, v = 2 • 10~ 9 . A force f n = F n (l + i)/u* n was used with F n real, 
nonzero only for n = 2,3, and chosen to give an input of energy but not 
helicity. Except where stated, the results shown are for the A = 26 simula- 
tion. Stationary time-averages were achieved by integrating over a period of 
more than 2900 large-eddy turnover times. In Fig. 1 are shown the energy 
spectrum and mean energy flux in the simulation. 
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n 

FIG. 1. Energy spectrum shows k~ 5 ^ 3 power law from 5th shell to 23rd shell. In 
the inset we show a constant mean energy flux in the inertial range. 

We first present evidence for the good decay of correlations of "spins". 
We define the spin-spin correlation functions Cxy(n\ m) = (X*Ym) — (X*)(Y m ) 
for X, Y = a, U. The results are shown in Figs. 2-3. It may be seen that the 
correlations decay quite rapidly in \m — n\, exponentially or as a large inverse 
power (> 7). This is an indication that the Gibbs measure of the hypothesis 
does not correspond to a critical point with a power-law scaling. For the 
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"1-dimensional" spin chain of the shell model a phase-transition would, in 
any case, require a long-range potential <1>, e.g. a pair interaction decaying 
by a small inverse power < 2 (e.g. see [10]). The hypothesis of a short-range 
potential therefore rules out any such critical behavior. 
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FIG. 2. C aa (n]m), Re CV[/( n ; m ) and Ira Cuu{ n \ m ) f° r n=12. 
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FIG. 3. Re C,ju(n;m) and Im C a u(n;m) for n=12. 

Further evidence against any critical behavior comes from a considera- 
tion of the thermodynamics. The Gibbs hypothesis implies that a "Gibbs 
free-energy" g should exist for the spin models, defined by a suitable "ther- 
modynamic limit". Thus, we may introduce "magnetic fields" p correspond- 
ing to the u-spins and h = h x + ih y corresponding to the ?7-spins, so that 
the concave free-energy is defined by 



g(P,h)= lim — \nZ n (p, h) 

n— >oo n 



(5) 



with the "partition function" 



Z n (p,h) : = (exp 



J2(pa k + Re(h*U k )) 



k=l 



>• 



(6) 



The absolute structure functions are proportional to the "partition func- 
tions" at h = 0: {\u n \ p ) = kn P ^ 3 Z n (p, 0). Thus, the existence of a thermody- 
namic limit, as implied by the Gibbs hypothesis, yields a power-law scaling of 
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the structure functions ~ k n with the anomalous exponent 5( p := Cp — | re- 
lated to the free-energy by g(p, 0) = 5£ p 4n A. However, the Gibbs hypothesis 
implies also a power-law scaling for the "phase structure-functions" Z n (0, h). 
Fig. 4 shows clean power-law ranges for these quantities, solid evidence in 
favor of the Gibbs hypothesis (but not a proof, since the thermodynamic 
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FIG. 4. Phase structure functions (a) Z n (0, 0, -0.75); (b) Z n (0, -0.5, 0). 

limit could exist even for a non-Gibbsian measure; e.g see [11]). In Fig. 5 we 
plot cross-sections of the Gibbs free energy, g(0, h x , 0) and g(0, 0,h y ). These 
appear to be smooth functions of their arguments. There is no evidence for 
any non-analyticity that would signal appearance of a phase transition. 
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FIG. 5. Gibbs free energy (a) g(0,h x , 0); (b) g(0,0,h y ). 

We now study the question of "translation-invariance" of the measures, 
by means of the single-site spin distributions. We plot in Figs. 6 and 7 
the distributions P(cr n ) and P(A n ), for different values of n in the inertial 
range. We see that these distributions collapse quite well, verifying the 
"translation-invariance" assumption. The first distribution is approximately 
exponential type P(cr) ~ (a/2) exp(— a\a\) with a ss 2.00 while the second 
is fit well by P(A) ss Cexp(-/3sin(A)) with C » 0.1, (3 w 1.2. The "infinite- 
temperature" model in [3] does not predict well either P(cr) or P(A). That 
approximation assumed a distribution P(A) uniform on the interval [— ir, 0] 
and yielded a multiplier distribution P(w) also compactly supported on a 
finite interval [u>-, w+]. However, the result in Fig. 6 (see inset) implies 
a distribution P(w) with two power-law regimes, ~ w 01 ^ 1 for w <^ 1 and 
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~ w' ^ 1 for w S> 1. This is inconsistent not only with [3] but with any 
independent spin model. Because of the power-law tail for large w, the 
moments which would give the anomalous scaling exponents for independent 
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FIG. 6. Distributions of a n , for n = 9 — 12. In the inset are the distributions of 
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FIG. 7. Distributions of A„, for n = 9 - 12. 

The distributions for the cr-spins appear symmetric under spin-flip a — > 
— (j around the K41 value a = 0. In fact, the exponential tails on the right 
and the left arise symmetrically from the same events in SABRA, which we 
call "defects". These are events in which the amplitude in a single shell, 
say, the nth, drops to a very low value, p n <^ 1. Intermittent bursts in the 
shell models are generally preceded by such events [8, 12], but defects do 
not need to appear in association with bursts. Since a n = ln(/9„/p n _i), the 
negative tail where a n — 1 comes from realizations with a "defect" in the 
nth shell, while the positive tail where a n 3> 1 comes from realizations with 
"defects" in the (n — l)st shell. The statistical distribution of such "defects" 
in a given shell n can be inferred from the constancy of P{u n ) at u n = 0. 
By a change to polar coordinates (p n ,@n), one finds P(p n ) ~ (const.)p n 
for p n 1. For events in the negative tail, p n <C 1 while p n -i ~ 1- In 
that case, a n ~ In p n and p n ~ e°™ . Therefore, by change of variables, 



d(p n ) 



P(Pn) 



(const.) 



2cr n 



for a n <§i — 1. An identical 
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argument for the positive tail using a n ss — ln(p„_i) and p n _i ~ e~ CTn gives 
-f (°"n) ~ (const.)e~ 2e7n for cr n S> 1. This argument assumes only that the 
amplitudes p n , p n ~i in the ratio are not strongly correlated, which could 
suppress the long tails. The distribution P(A) also exhibits a symmetry 
A — > 7T — A around the K41 solution with A = —tt/2. This symmetry can 
be expressed as U — > —U* and is seen as well in the vanishing of Im Cjju an d 
Re C a u in Figs. 2-3 and in the near symmetry of g(0, h x , 0) in Fig. 5. In fact, 
this is an exact symmetry of the dynamics, broken only by our forcing. If 
u n is a solution of SABRA with a force f n , then —it* is a solution with force 
— /* and under this transformation U n — > — U* for all n. This symmetry 
should be restored for large n, similar to restoration of isotropy in 3D. 

The distribution P(a) cannot be exactly symmetrical under the spin-flip 
a — > —a. There must be a non-vanishing mean (a) or a "magnetization" , due 
to the fact that Kolmogorov 1941 mean-field scaling of pth-order structure 
functions is not exact at p = [13]. Indeed, by the Gibbs hypothesis, 
the magnetization can be obtained from the thermodynamic formula (a) = 

— If -In Fig. 8 we plot (a n ) vs. n for the two simulations of SABRA 
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FIG. 9. (Re(U n )) vs. n with — dg/dh x \ p h=0 and of (Im{U n )) vs. n with — dg/dh y \ p h=Q . 

with N = 22 and N = 26. We see that there is a slight breaking of trans- 
lation symmetry. This is likely due to finite-size effects and becomes much 
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smaller for the N = 26 shell simulation than for N = 22. In the same 
figure is plotted a straight horizontal line for the prediction (a) ~ —0.0757 
obtained from the derivative of g. We see that the agreement is quite sat- 
isfactory. It is also apparent from Fig. 7 that a non-zero expectation (U) 
occurs, non-invariant under conjugation U — > U* . In fact, a negative value 
(sin(A)} < is required in SABRA for a forward cascade of energy [3, 4]. 
We plot in Fig. 9 the "magnetizations" (U n ) vs. wavenumber n. Again 
we see a finite-size breaking of translation-symmetry, which lessens going 
from N = 22 to N = 26. The numerical values of these "magnetizations" 
are also predicted by derivatives of the Gibbs free energy, giving a large ex- 
pectation (sin A) « —0.489 but a small value (cos A) = 0.00814, consistent 
with approximate symmetry U — > —U*. These results are plotted in Figs. 8 
and 9 as horizontal lines and obviously give satisfactory agreement with 
the direct measurements of magnetizations. The results indicate that there 
is "spin-ordering" in the turbulent systems, which breaks discrete spin-flip 
symmetries. Peierls-type arguments [14], including rigorous versions such 
as Pirogov-Sinai theory [15], indicate that discrete symmetries cannot be 
spontaneously broken in "1-dimensional lattice" systems such as the shell 
model, if the interaction potential is short-ranged. A 1-dimensional Gibbs 
distribution with non-zero magnetization for a symmetric, short-ranged po- 
tential is unstable to formation of domain walls. We expect the spin order 
here to be not "ferromagnetic" but "paramagnetic", arising from explicit 
symmetry-breaking terms in the potential <& of the Gibbs measure. 

4 Discussion and Conclusions 

In this work, we have shown that Kolmogorov multipliers and dynamical 
phases in the shell model must be correlated from shell-to-shell. We have also 
presented evidence that these "spin" variables, while not independent, are 
distributed according to a translation-invariant Gibbs measure. Of course, 
this has certainly not been proved in this work. There are non-Gibbsian 
measures on one-dimensional lattices, such as that of Schonmann [16], which 
have both exponential decay of correlations and a thermodynamic limit of 
the Gibbs free energy (which, however, is there non-analytic). On the other 
hand, that type of example has been shown to be Gibbs also in a somewhat 
generalized sense [17]. The hypothesis that the distribution of suitably de- 
fined "spins" in the shell model is Gibbsian with a summable potential has 
testable consequences, some of which we have verified in this work. 

If the hypothesis is true, then it is possible to recover the potentials from 
the finite-shell marginal distributions. Indeed, if P n ,...,n+N-i(£n, ■•■>£n+iv-i) 
is the probability density of the "spins" at shells n,...,n + N — l in the 
inertial range, then 

N-l 

lnP ni ... ;n+iV _i(£ n ,...,£ n+J v--i) = - $n+k(Z) + °( N ), (7) 

fc=0 
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whenever the potential is absolutely summable. See Proposition 2.46 in [18]. 
More directly, the potentials may be derived from conditional probabilities of 
the spins via the Mobius inversion formula [19]. For example, if the "spins" 
are distributed by a Markov chain or nearest-neighbor Gibbs measure, then, 
up to constants, the 1-body interaction is 

4% n ) = -mT^^lf ) -mT n+1 | n (C|£n) (8) 
and the 2-body interaction is 

<Ll(£n,£n-l) = -lnT^^l^x) 

+ lnT n | n _ 1 (£„|0 + lnT n | n _ 1 (C|£n-i). (9) 

Here T n | n _ 1 (^ n |^ n _i) is the transition probability of the Markov chain and 
is a constant reference value of the spin. As this example makes clear, 
the potentials of the Gibbs measure can be recovered, in principle, from 
conditional probabilities obtained in numerical simulations of the dynamics. 

If true, it is intellectually interesting that probability measures arising 
from turbulent dynamics may be Gibbsian, in a suitable spin representation. 
However, more importantly, it gives one some new tools that one may apply 
to the turbulence problem. For example, it offers a new route to calculate 
the scaling exponents as a "free energy" g(p, h) of a one-dimensional spin 
system. Furthermore, formulas for the conditional probabilities of small- 
scale modes given the large-scale ones provided by the Gibbs hypothesis 
may be very useful in carrying out large-eddy simulations of turbulence 
[20] . Some preliminary tests of these ideas have already been carried out for 
the shell models and will be reported elsewhere [6]. 

Similar results as discussed here hold for Navier-Stokes dynamics [21], 
using the velocity increments advocated by Kolmogorov [1] . Multipliers and 
spins may also be defined by a representation of the Navier-Stokes equa- 
tions in terms of orthogonal wavelet bases [22]. In that formulation, the 
dynamics resembles the shell model on the dyadic Cayley tree [23]. Gener- 
alizing [3] , the invariant measure of the latter should be a Gibbs measure on 
a Bethe lattice ([24], Ch.4.) Whereas standard shell models corespond to 
"1-dimensional" spin systems, the spin systems on the Bethe lattice are effec- 
tively "infinite-dimensional". Nevertheless, the statistics in the shell model 
and in Navier-Stokes should be qualitatively similar, as present evidence 
suggests that the Gibbs distributions for both are in the high-temperature, 
unique-phase regime. 

A Appendices 

A.l Stabilization of the K41 Solution 

We show in this first appendix that the K41 solution would be stable if 
multipliers were independent for distinct shells (or, in fact, even if they 
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were independent for shells a finite distance apart). We show that this 
follows from a set of exact dynamical constraints on the multiplier variables 

£,n = (w n ,A n ). 

For this purpose, we must transform the equations of motion into those 
variables. In terms of p n ,9 n , the SABRA dynamics becomes 



,2/3 
Pn Kn 

Pn A 



and 

,2/3 

On ~ — 



Pn+2Pn+l . a , ,Pn+lPn-l ■ A . Pn~lPn-2 . A 

a sm A n+2 + b sin A n+ i + c sin A n 

Pn Pn Pn 

. kn Re(/ n e l6n ) 2 / in \ 

+ »K (10) 

Pn 



Pn+2Pn+l * , ,Pn+lPn-l A Pn-lPn-2 A 

a cos A n+ 2 + 6 cos A n+ i — c cos A n 

Pn Pn Pn 



Pn 

Going over to the scale- local variables w n , A n , this becomes 
w n = w n [U n (w,A) - U n -i(w, A)] := W n (w, A) 

and 

A„ = -V^w, A) + V n -i{w, A) + V n - 2 (w, A) := Z n (w, A) 

where U n ,V n are the righthand sides of (10), (11), respectively, expressed 
in terms of w n ,A n . Thus, considering just the inertial-range part of the 
dynamics, 

k 2/3 



U n (w,A) = • [a • w n+2 smA n+2 ■ w 2 n+l w n w n -i 



n-2 



+ b ■ w n+ i sin A n+ i • w n -i + cw n 1 sin A n ] w fc (12) 



fc=i 



and 



fe 2/3 

T4(U), A) = • [a • W n+2 COS A„ +2 • U^+iWnWn-l 



n-2 



+ 6 • w n+ i cos A n+ i • w n -i - cw n 1 cos A n ] w k (13) 

fc=i 

In terms of these variables, the dynamics appears highly nonlocal in scale, 
because of the product IIfc=i w k- 

If P n {S,n) is the distribution of £ n = (w n ,A n ), then it is straightforward 
to show that 

BtPniin) = [<W„|e„>P„(C„)] " [<2„|&,)P„(£„)] , (14) 
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where (W n |£ n ), {Z n \£ n ) are conditional averages for fixed £ n . These equa- 
tions are exact, but not closed in terms of P n . However, if we assume that 
the statistics of the model are given by a product measure Y\ n Pn(£,n) — and 
with only that assumption — then we obtain closed equations for the -P n 's. 
In the inertial range, these equations are of the form (14) with 

^2/3 ra-4 

(W n \Cn) = Y[{w k ) [{A n w 3 n - B n w 2 n + C n w n ) + (D n w 2 n - E n ) sin A n ] 
fc=l 

k 2/3 n " 4 

(Z n \£ n ) = Y[{w k ) [{F n w 2 n + G n w n + H n ) - {K n w~ l + L n w n ) cos A„] 
k=l 



where 



B n = -a(w n+2 S n+2 )(w n+1 )(w n -i)(w n - 2 )(w n - 3 ) 

c S _i 

C n = 6(w n+ iS , n+ i}(w n _i)(u) n _ 2 }(w n _3) - TTn(— )( w n-3) 

E n = -c(w n - 2 }{Wn-3) 
F n = -^J^{w n +lC n+1 )(w n - 1 ){w n -2){w n -z) 

G n = -a{w n+2 C n+2 ){w 2 l+1 ){w n -i){w n - 2 ){w n -r i ) 

c C n -i 

H n = -b(w n+1 C n+1 ){w n - 1 ){w n - 2 ){w n -2,) ~ \h-K(— — )(w n -3) 

+^K-iC„_ 1 )K_ 3 ) - ^s(^) 
K n = -c(w n _ 2 )(w n - 3 ) 

In these expressions S n = sinA n , C n = cosA„ for all n. Note that the 
resulting closed equations for the P n 's are nonlinear integro-partial differen- 
tial equations, since the averages in the above expressions are over the P n 's 
themselves. 

If we assume further that all the "spins" £ n are identically distributed, 
i.e. P n = P for all n, then the above equations for P n for each n reduce 
to the same equation for P, after a change in the time-scale by a factor 

4rM n_4: 

9tP(o = ~ mop(0] - ^ mono] , (15) 
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with now 

W(£) = (Aw 3 - Bw 2 + Cw) + (Dw 2 - E) sin A 
Z(£) = (Fw 2 + Gw + H)- (Kw~ l + Lw) cos A 

where 

B = -a(wS)(w 2 )(wf 
C = b(wS)(wf-^£)(w) 

2 



E = -c{w) 



G = -a{wC)(w 2 )(wf 
H = -b( W C)( W f - ^£)M + ^(wC)(w) - 

K = -c{w} 2 

L = "ipH ! - j^<»*>M a - 

Using a + 6 + c = 0, it is easy to verify that P±(£) = 6(w — 1)6 (A. =F §) are 
exact time-independent solutions of (15). The presence of these solutions 
is due to the well-known existence of exact, steady-state "K41" solutions of 
the SABRA model, of the form = ±iAkn X ^ for any choice of the real 
constant A > 0. The solution has a backward energy transfer to low 
wavenumbers, while u~ has forward transfer to high wavenumbers. 

It is interesting that P-(£) is linearly stable under the dynamics (15), 
while P+(£) is unstable. In fact, the linearization of (15) is 

9tSP(0 = -^[w ± (0sp(0]-^ k [z±(0sp(0] 

with 

W±(0 = (A±w 3 - B±w 2 + C±w) + (D±w 2 - E±) sin A (17) 
Z±(£) = ~(K ± w" 1 + L±w) cos A (18) 

where 
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F± = G± = H± = 0,K± = -c,L± = -[— + 



A 2/3 A 4/3 

and 5W(£),8Z(£) have the same form as W(£), Z(£) but with coefficients 
5 A, 8B, etc. that can be obtained by linearizing the corresponding coeffi- 
cients A, B, etc. Now the essential fact is that £± = (1, ±|) is an unsta- 
ble/stable fixed point of the dynamical system (w, A) = (W±(£), Z±(£)), for 
± respectively. This is easily verified directly from the equations (17), (18). 
For example, the linearization about those fixed points is 





dZ- 


dw 


dA 


dW- 


dW~ 


dw 


dA 


d~Z^ 


az+ 


dw 


dA 


dW+ 


dW+ 


dw 


dA 



(l-^Ts) (a-c) 
/(2a + 6) + 2|+p 

V i}-w*) a+ {}-w-z 



When a > and 6, c < 0, then we see that both eigenvalues are positive 
for the + fixed point, and the second eigenvalue is negative for the — fixed 
point. The first eigenvalue of the — fixed point is also negative when 

2a + 56 > -(2a + 6)A 2/3 , 

which imposes a condition on the coefficients. For example, with the com- 
mon parameterization a = 1,6 = — e, c = e — 1, a region in the (e, A)-plane 
is selected specified by 2 — 5e > (e — 2) A 2 / 3 . Along the curve A = 1/(1 — e) 
for which the second invariant is helicity-like, the condition is that A should 
lie in an interval (1,A*) with A* ~ 2.467. In particular, the standard case 
A = 2, a = 1,6 = c = —1/2 which we simulated in this paper lies in this 
region. For all parameter values satisfying the above condition £ + is linearly 
unstable, and linearly unstable. 

Let then Bs (£±) be a disk of radius 5 centered at £± in the strip (0, 00) x 
[— 7r,7r] of the (w, A)-plane, and let Bg(^±) be its complement. Since the 
integral over B^((,±) of the second set of terms in (16) vanishes identically, it 
follows that the perturbation term j B c^ ± ^ d£, 5P(£) satisfies the same equa- 
tion as it would under the flow of the vector field in (18)-(17). In that case, 
for the + sign, J B c^ + \ d£ $P(0 increases in time, whereas for the — sign 
Jgc(£ ) dl; 5P(£) decreases in time, for any 5 > 0. We therefore see that 
solution -P+(£) is linearly unstable, whereas the solution P-(£) is linearly 
stable. In fact, a direct integration of the full nonlinear equation (15) shows 
that P-(£) is the global attractor for all initial probability density functions 
P(0- See [6]. 

This argument, given here assuming perfect independence of "spins" at 
different shellnumbers, can be generalized assuming a correlation of finite 
range r. In that case, one can develop a similar equation for the r-spin 
distribution P r (£ n , £ n +i> ■■■■> £n+r)- It is found in the same manner that a 
delta-function at the K41 fixed point value £_ in all shells is the unique, 
global attracting solution. For details, see [6]. 
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A. 2 Markov Chain Models of Multipliers 

In this second appendix we briefly discuss Markov chain models for the 
multipliers. We shall assume that £ n = (w n , A n ) has statistics derived from 
a stationary Markov chain with single-shell distribution P(£ n ) an d forward 
transition probability T(£ n+ i|£„), n = 1,2,3, ... Because of the stationarity 
assumption, these functions do not depend upon shellnumber n. We first 
discuss the relation between the structure-function scaling exponents and 
the eigenvalues of certain "transfer matrices" . This is a particular example 
of the connection discussed in the text between scaling exponents and free 
energy functions. We next discuss the scaling properties of the mean fluxes 
of conserved quantities. In particular, we show how Markov chain models 
can yield a joint cascade of both energy and helicity. 

To calculate the absolute structure functions, it is easiest to consider the 
time-reversed Markov chain, with backward transition probability 

~ _ T(Z n+1 \£ n )P(Z n ) 

1 {Qn\Qn+l) — ^TZ ; • 

P{U+i) 

Then, using p p n = Uk=i < = 112=1 e p<Tfe , 

{&) = J <%n+i P(Z n+ i) J d^ n e^T(£ n |£ n+1 ) ■ • • J d£i e pCTl f (616) 

= j di n+1 P(Ca + l) J d£l 7f p) (6|£n+l) (19) 

where the "transfer matrix" TL) is defined by 

T (p )(Zn\Z n+ i):=<? ffn T(S n \S n+1 ) 

and is its n-fold convolution or matrix product. For large n, we have 
asymptotically that 

~ ^ p) R( P )^l)L {p) (Cn + l) (20) 

where is the principal eigenvalue of T( p ) and R^,L^ are correspond- 
ing right and left eigenfunctions. By the Perron-Frobenius theorem, these 
eigenvalues are real and non-negative, and the left and right eigenfunctions 
may also be chosen to be non-negative. Thus, we find that 

for n — > oo, where {Lr p \) = J d£ P(£)L( p )(£) and Rr p \ = j d£ Rr p \(£). In this 
way, we obtain the relationship between the scaling exponents of structure 
functions and eigenvalues of the transfer matrices as 

Cp = | - log A (/i(p)). 
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Exact moment constraints require that /ir^ = 1, so that ^3 = 1 within a 
Markov chain model. For details, see [6]. 

We now consider the asymptotics of mean fluxes of the inviscid invari- 
ants. In the SABRA model, the energy flux is given by 



H71 — O^n'^ra+l + ck n —iT n 



and the helicity flux by 



u { (X \ n 

= -2a - (k n T n+1 - K-xTn). 



. c 

Here T n is a triple velocity product 

T n = Im(u* +1 n„'U n _i) = -^-p n+ ip n p n _isinA„ + i 

tin 

In terms of the scaledocal variables w n , A n , n = 1,2, 3, ... 



n-l 



1 

Ilf = - [-a • w n+2 sin A n+2 • • + c • w n+ i sin A n+i • w r 2 J J| w| 

fc=i 

and 

nf = — (-) [w n+ 2 sin A„ +2 • wl +1 ■ wl - w n+ i sin A n+ i • u£] 

A;=l 

We now evaluate the mean values of these flux variables. If ^2(^^+25 £n+i) — 
T{£ n +i\£n+2)P{in+2) is the joint distribution of "spins" £„ +2 ,£n+i at two 
successive shells, then 

<nf ) = i 

-a y d£ n+ 2 y (i^ n+ l P 2 (Cn+2,Cn+l) • ^n+2SinA n+ 2 • • 



+c y d£ n+ i y d£ n P 2 (£ n +i,£n) • Wn+i sin A n+ i -u^ y d£i T™ } H^llCn) 



(21) 



and 

2a /a\" 



y d£n+2 y P2(£n+2,£n+l) ■ W„+ 2 sin A n+2 • U>^ +1 ■ y d£i Tj^ 

y <% n P 2 (tn+i,tn) ■ w n +i sin A n+1 ■ wl J d£i T™ } 1 {Zi\£ n )\ (22) 

To calculate the mean energy flux asymptotically for large n, it suffices 
to use the previous asymptotic expansion (20) for p = 3. This gives 

( n n) ~ -j(W(3) ~ c)^- ) 1 (w 2 smA 2 ■ wlL {3) { Wl , Ai))^ (3) 
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as n —y oo. The only way that the energy flux can be asymptotically con- 
stant, i.e. independent of shellnumber n, is if ^(3) = 1. This is another 
argument for that constraint, independent of that given in [6]. In that case, 
we obtain finally that 



( n n) j(a - c)(w 2 smA 2 ■ wjL^iwx, A^R^ 

asymptotically for n — ► 00. 

However, for helicity flux, the contribution from the leading-order term 
in the asymptotic expansion (20) gives zero identically, because 



w 



l+2 sin A n+2 • uv 2 t+1 L (3 ) (w n+1 , A n+ i)} = (w n+1 sin A n+1 ■ ^L (3) (w n , A n )}, 



by stationarity of the Markov chain. Thus, a non-vanishing contribution is 
obtained only from the next order in the asymptotic expansion for large n, 

T (3)(£n\£n+l) ~ ^3)^(3) (Cl)^(3)(Cn+l) + ^(3)^(3) (&) L '(3) (Cn+l) 

where // 3 n is the subleading eigenvalue of T( 3 ) (i.e. the complex eigenvalue 
with next largest magnitude |/^ 3 )| after |/U( 3 )|), and R'^y L'^ the correspond- 
ing right and left eigenfunctions. In that case, 

This flux is constant precisely when //^ = c/a. In that case, 



(nf> - y c (a-c)(w 2S mA 2 -w 2 1 L[ 3) (w 1 ,A 1 ))R' (;i) 

asymptotically for n — > oo. 

When c/a < 0, a non-zero helicity flux is ruled out in an independent 
multiplier model by the realizability inequality (u? 3 ) > 0. However, in a 
Markov chain model, there is no such constraint, because the subleading 
eigenvalue fi'^ may easily be negative. As a concrete example, consider the 
Markov chain with single-shell distribution 

P{tn)=Ce-^ A -P W { Wn ), 

where Pw{w n ) is any density on the interval (0, oo), and with transition 
probability density 

T(Cn\Cn+i) = P(£n) ■ 1 + (^) sgn(cos A n )sgn(cos A n+i ) 

The latter kernel is non- negative when |c/a| < 1. In fact, it is a rank 2 
operator and can be written as 

f (£„|£„+i) = i?(£ n )L(£„ + i) + (-) R'(Ca)L'(C n+ i) 
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where 

R(tn) = P(£ n ), L(£ n ) = 1, R'((; n ) = P(£„)sgn(cos A n ), L'(Cn) = sgn(cos A n ). 

It is easily checked that these four vectors are a bi-orthonormal set, and can 
be completed to a bi-orthonormal basis. It follows directly that 

j d£ n ?(£ n |£ n+1 ) = 1, 

so that T is a transition probability density, as claimed. Also P is its 
invariant distribution, because 

j da n+1 r(£ n |£ n+1 )P(£ n+1 ) = 

Furthermore, if J dw w 3 Pw(w) = 1, then by construction the only two non- 
vanishing eigenvalues of T( 3 )(£ n |£ n+ i) = wlT(£ n \£ n+1 ) are // (3) = 1 and 
^(3) = c / a ' w ith 

#(3)(£n) = U^&O, L (3) (e„) = l 
#(3)(&0 = w n P (£n) s g n ( cosA r t )i L '(3){€n) = Sgn(cOsA n ). 

This example is not completely realistic as a statistical model for the shell 
dynamics, because the amplitude multipliers w n , n = 1,2, 3, ... form an i.i.d. 
sequence. Hence, it is not consistent with the long power-law tail ~ w~ 3 
which we have observed in the numerical simulation. It is also does not 
have a non- vanishing flux of helicity because, unfortunately, the expectations 
(102 sin A2 ■w\L'^ j {w\, Ai)) = R{z) = in this model. However, it illustrates 
how it is possible to get the subleading eigenvalue c/a < within a realizable 
Markov chain model. 
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